Banks of templates for directed searches of gravitational waves 
from spinning nuetron stars 
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We construct efficient banks of templates suitable for directed searches of almost monochromatic 
gravitational waves originating from spinning nuetron stars in our Galaxy in data being collected by 
currently operating interferometric detectors. We thus assume that the position of the gravitational- 
wave source in the sky is known, but we do not assume that the wave's frequency and its derivatives 
are a priori known. In the construction we employ simplified model of the signal with constant 
amplitude and phase which is a polynomial function of time. All our template banks enable usage 
of the fast Fourier transform algorithm in the computation of the maximum-likelihood ^-statistic for 
nodes of the grids defining the bank. We study and employ the dependence of the grid's construction 
on the choice of the position of the observational interval with respect to the origin of time axis. 
We also study the usage of the fast Fourier transform algorithms with non-standard frequency 
resolutions achieved by zero padding or folding the data. In the case of the gravitational-wave 
signal with one spindown parameter included we have found grids with covering thicknesses which 
are only 0.1%-16% larger than the thickness of the optimal two-dimensional hexagonal covering. 

PACS numbers: 95.55.Ym, 04.80.Nn, 95.75.Pq, 97.60.Gb 



I. INTRODUCTION AND SUMMARY 

Rotating neutron stars in our Galaxy are expected 
sources of almost monochromatic gravitational waves 
which arc looked for in the data being collected by cur- 
rently operating ground-based interferometric detectors 
LIGO Q, Virgo [3, GEO600 j|, and TAMA 0]. In 
the present paper we consider the problem of construc- 
tion of efficient banks of templates needed to detect 
gravitational- wave signals originating from spinning neu- 
tron stars. We assume that the problem of detection of 
the signal and of estimation of its parameters is based 
on the maximum-likelihood (ML) principle and we also 
assume that the noise in the detector is Gaussian and 
stationary (the detailed exposition of the ML detection 
in Gaussian noise can be found e.g. in Chapter 6 of the 
monograph see also [f| and the review article Q). 
Data analysis tools and algorithms needed to perform, 
within the ML approach, an all-sky search for almost 
monochromatic gravitational- wave signals (i.e. the search 
which assumes that the position of the source in the sky 
is not known) were developed in detail in the series of 
papers [§-[T2j (see also Refs. 0,0)- in our paper we re- 
strict ourselves to directed searches for almost monochro- 
matic signals, i.e. we assume that the position of the 
source in the sky is known. We do not assume however 
that the frequency and the spindown parameters of the 
gravitational-wave signal are known. Directed searches 
should thus be distinguished from targeted searches, in 
which it is additionally assumed that the frequency and 
the spindown parameters of the gravitational wave are 



also known (see Refs. [HI, [l6| for discussions of statistics 
by means of which one can test whether data contains 
such signal). Several targeted searches were already per- 
formed with data collected by the LIGO and GEO600 
detectors G3-II3. 
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In the ML approach one considers the likelihood ratio 
A[x; 9} which is a function of the data x and the param- 
eters of the gravitational-wave signal we are looking 
for. Detection of the signal relies on the computation of 
A[x; 0] maximized over all possible values of the param- 
eters 6 and comparing this maximum with a threshold. 
In the case of directed searches the unknown parameters 
6 can be divided into two groups, 6 = (A,£). The first 
group A consists of four extrinsic or amplitude param- 
eters: an overall amplitude of the waveform, its initial 
phase, the polarization angle of the wave, and the incli- 
nation angle of the star's rotation axis with respect to 
the line of sight. The second group £ contains intrin- 
sic or phase parameters: the frequency of the wave and 
the spindown parameters. Maximization of the A with 
respect to amplitude parameters A can be done analyt- 
ically (by solving the set of ML equations dA/dA = 
with respect to A) and the J- ' -statistic is defined as the 
logarithm of the likelihood ratio A after replacing in A 
the amplitude parameters A by their ML estimators A: 
T\x;£\ := lnA[x; A, £]. Maximization of the ^-statistic 
over the phase parameters £ can be done only numeri- 
cally. To find this maximum one needs to construct a 
bank of templates in the space of the parameters £ on 
which the ^-statistic depends. The bank of templates is 
defined by a discrete set of points, i.e. a grid in the pa- 
rameter space chosen in such a way, that for any possible 
signal there exists a grid point such that the expectation 
value of the ^-statistic (which is a random variable as 
it depends on the detector's noise which is a stochastic 
process) computed for the parameters of this grid point 
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is not less than a certain fixed minimal value (assuming 
that the minimal value of the signal-to-nise ratio is also 
a priori fixed). 

In the series of papers [814121] it was argued and checked 
by numerical simulations that in the case of an all-sky 
search in the construction of the bank of templates one 
can employ a simplified model of the gravitational-wave 
signal, the so called linear phase model, in which the am- 
plitude of the signal is assumed to be constant and the 
signal's phase is a linear function of the unknown param- 
eters (this model was introduced in Sec. V B of Ref. [9]). 
In Sec. V E of [9J it was checked that the linear model 
reproduces well the covariance matrix (defined as the in- 
verse of the Fisher matrix) for the ML estimators of the 
signal's parameters of the exact gravitational-wave sig- 
nal. Also in Ref. @ (Sec. V D and Appendix C there) 
the polynomial phase model was introduced in which the 
signal's amplitude is constant and the phase is a polyno- 
mial function of time. It was found (in Sec. V E of [§]) 
that the polynomial phase model reproduces very well the 
covariance matrix of the signal's parameters of the exact 
model in the case of directed searches. This indicates 
that the polynomial model can be used in the construc- 
tion of the banks of templates for directed searches and 
this model is accepted in the present paper. 

In the current paper we assume that the observational 
interval is of the form (ti—T /2; t[+T /2) (where T is the 
length of observation time), and we study the dependence 
of the construction of banks of templates on the choice of 
the parameter ti (or its dimensionless version \i U/T ) 
which fixes the position of the observational interval with 
respect to the origin of time axis. 

The organization and the main results achieved in the 
paper are as follows. In Sec. II we introduce the poly- 
nomial phase model of the gravitational-wave signal. We 
consider here the phase with only one spindown parame- 
ter included. For this model we compute the ^"-statistic 
and its expectation value in the case when the data con- 
tains the gravitational-wave signal. This expectation 
value depends [see the crucial Eq. (I2.33[) below] on the 
signal-to-noise ratio p and on the value of the autocovari- 
ance function Cq ) of the ^-statistic (the subscript 
'0' indicates that the autocovariance is calculated in the 
case when data is a pure noise) computed for the intrin- 
sic parameters of the template (£) and the gravitational- 
wave signal (£'), respectively. The signal-to- noise ratio p 
we can not control, therefore to construct bank of tem- 
plates one needs to choose some minimum value C min of 
the autocovariance function Co and look for such a grid 
of points that for any point in the intrinsic parameter 
space there exists a grid node £ such that the autoco- 
variance Co (£,£') computed for the parameters £ and 

is not less than C m ; n . The autocovariance Co(r) (for 
polynomial phase model it depends on £, only through 
the difference t := £ — £') can be expressed in terms of 
Fresnel integrals or one can use an approximate formula 
for Co by taking the Taylor expansion (up to the second- 
order terms) of Cq around its maximum at t = 0. In Fig. 



Q] and Table Q] we compare these two ways of computing 
the autocovariance. 

In Sec. Ill we consider banks of templates neces- 
sary to perform detection of almost monochromatic 
gravitational-wave signals with polynomial phase. We 
first formulate the problem of constructing bank of tem- 
plates as a problem of finding optimal covering of the 
signal's parameter space by means of identical ellipses 
(defined as isoheights of the autocovariance function Co 
of the ^-statistic) and introduce some mathematical no- 
tions related with coverings. We are interested in such 
searches for almost periodic gravitational- wave signals for 
which the number of grid points in the parameter space is 
very large and the time needed to compute the ^-statistic 
for all grid nodes is long. Then it is crucial to use in the 
computation the fast numerical algorithms. Because the 
computation of the ^-statistic involves calculation of the 
Fourier transform, one would like to use the fast Fourier 
transform (FFT) algorithm. As is known, the FFT algo- 
rithm computes the values of the discrete Fourier trans- 
form (DFT) of a time series for a certain set of discrete 
frequencies called the Fourier frequencies. Thus it will 
be possible to use the FFT algorithm in computation of 
the ^-statistic, provided the grid points will be arranged 
in such a way, that the frequency coordinates of these 
points will all coincide with the Fourier frequencies. We 
have constructed two different families of grids which ful- 
fill this requirement. Our constructions were motivated 
by a simple observation (see Fig.[T} that the shape of the 
autocovariance ellipse strongly depends on the value of 
the parameter Xi- the larger is, the more along ated 
(along the frequency axis) this ellipse is. The usage of 
alongated ellipses suggests that the ^-statistic could be 
computed for smaller number of frequency values than 
the standard DFT algorithm computes. Therefore we 
have employed some modifications of the DFT algorithm 
leading to non-standard frequency resolutions achieved 
by zero padding or folding the data (this is discussed in 
Appendix A). Details of the grids construction are pre- 
sented in Appendix [Bj where also the components of the 
basis vectors spanning the grids are given (in Tables IIIH - 

E}. 

The construction of efficient banks of templates for 
gravitational- wave searches was recently discussed in Ref. 
[22j | (see also [23h26| ). where random template banks 
and relaxed lattice coverings were considered. As ex- 
plained above we are interested in searches involving data 
streams so long, that the time performance of the search 
crucially depends on the ability of using the FFT al- 
gorithm. This enforces the above-mentioned constraint 
which is not fulfilled by the grids considered in Ref. [22| ■ 
Therefore our work can be considered as being comple- 
mentary to the studies performed in Ref. |22j. The grid 
fulfilling the constraint was constructed in Sec. IV of Ref. 
[l2| in the case of all-sky search for almost monochro- 
matic gravitational waves. In the search considered in 
the signal's parameter space is 4-dimensional and the 
polynomial phase model can not be employed. However, 
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in Sec. IV C of [12(, as an example of application of gen- 
eral algorithm devised in [l2| to construct constrained 
grids, the two-dimensional grid for searches of signals 
with polynomial phase is considered. This grid has thick- 
ness equal to ~1.8, whereas grids constructed in the cur- 
rent paper have thicknesses of ~1.2 (all grids found by 
us have covering thicknesses which are only 0.1%-16% 
larger than the thickness of the optimal two-dimensional 
hexagonal covering). 

In Sec. IV of the current paper we present the results of 
numerical simulations we performed to study the quality 
of grids constructed in Appendix B. These results are 
contained in Table |TTJ From our simulations it follows 
that when considering the efficiency of a grid it is not 
enough to worry only about the thickness of the covering 
related with the grid. We see from Table [TT| that for 
lattices which have the same thickness of covering, but 
one of them has basis vectors which are orthogonal to 
each other, the time used to compute the ^-statistics is 
at least two times shorter than the corresponding time 
for the lattice with non-orthogonal basis vectors. 



II. AUTOCOVARIANCE FUNCTION OF THE 
^-STATISTIC 

We assume that the noise n in the detector is an ad- 
ditive, stationary, Gaussian, and zero-mean continuous 
stochastic process. Then the logarithm of the likelihood 
function is given by 



In K[x] = (x\h) - 1 (h\h), 



(2.1) 



where x denotes the data from the detector, h is the 
deterministic signal we are looking for in the data, and 
(• | • ) is the scalar product between waveforms defined by 



{hx\h 2 ) :=4Re 



hi{!)h* 2 (f) 
Sn(f) 



(2.2) 



where " stands for the Fourier transform, * denotes com- 
plex conjugation, and S n is the one-sided spectral density 
(defined for frequencies < / < +oo) of the detector's 
noise n. 

We are interested in almost monochromatic signals, i.e. 
such signals for which the modulus of the Fourier trans- 
form is concentrated (for frequencies / > 0) around some 
'central' frequency f c > and S n is a slowly changing 
function of / in the vicinity of the frequency f c . If both 
waveforms h± and hi in Eq. (|2.2[) have Fourier transforms 
concentrated around the same frequency f c , then we can 
replace S n (f) in the integrand of (|2.2[) by S n (f c ) and, 
after employing the Parseval's theorem, approximate the 
scalar product by 



(/ii|/i 2 ) = 



ti+T /2 

' hi(t)h 2 (t)dt=-^^(h 1 h 2 ) 



Sn(fc) 



ti-T /2 



Sn(fc) 



Here (ti — T D /2; ti + T Q /2) denotes observational interval, 
so T Q is the length of observation time, and tj — T a /2 is 
the moment at which the observation begins. The time 
averaging operator ( • ) is defined by 



U+To/2 



(h) := 



T 

- 1 o 



h(t) dt. 



(2.4) 



*i-T /2 



Using the formula (|2.3I) we can write the log likelihood 
ratio from Eq. (|2.ip as 

lnA[x] -iSh( {xh) ~l {h2) )- (2 - 5) 

We restrict to almost monochromatic signals with fre- 
quency drift given by the linear-in-time relation, 



/(t) = /o + M 



(2.6) 



so fo and /o is the frequency and the first derivative of the 
frequency with respect to time, respectively, both taken 
at the time t = 0. The instantaneous frequency is related 
to the phase \& of the signal by 



/(*) 



1 dtf(f) 
2tt dt 



(2.7) 



Making use of Eqs. (|2.6[) and (|2.7[) . we obtain the time 
dependence of the phase \I>: 

*(t) = 27r^ t /(i')d<' + *o = 27r (f 1 + I f t 2 ^j + 

(2.8) 

where is the value of the phase $ at the time t = 0. 
It is convenient to define the dimensionless parameters 

oj :=2TTf T , Wi :=tt/ T o 2 . (2.9) 

Using these parameters we can write the phase ^ in the 
form 



where 



*(£;*o,wo,wi) = $(t;wo,Wi) + * , 



t ft 
<P{t; w , U\) = w — + wi I — 



(2.10) 



(2.11) 



We further assume that the gravitational-wave signal h 
we are looking for has a constant amplitude ho , so it can 
be written in the form 

h(t; h , * , wo, wi) = ho sin ($(*; w , wi) + *o) ■ (2.12) 

Let us collect the parameters of the phase $ into a 2- 
dimensional vector £, 



(2.3) 



(2.13) 
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Then the gravitational-wave signal h, Eq. (|2.12[) . we can can represent the time average (h 2 ) as 
shortly write as 



h{t; h , o, £) = h sin ($(t; £) + tf „) 



(ft 2 ) = i(/i 2 + hi) + /ii/i 2 (sin2$) + ~(/jf, - /i 2 J(cos2$). 

(2.14) ( 2 - 17 ) 
For observation times longer than few hours and for al- 
most monochromatic signals with frequency of the order 
of hundreds or thousands hertz, we can make approxima- 
It is easy to maximize the likelihood ratio (|2.5p for tion 



the signal (|2.14l) with respect to the parameters ho and 
^o- To do this it is convenient to introduce the new 
parameters hi and h 2 : 



hi := ho cos^Oi h 2 :— ho sin^o, (2-15) 
and rewrite the signal h in the form 

h{t;hi,h 2 ,£) = hism$(t;£) + h 2 cos$(t;£). (2.16) 



(sin2$)=0, (cos2$)=0. (2.18) 
Then Eq. (12.171) can be approximated by 

{h 2 ) = \{h\ + hl). (2.19) 

With the aid of the equality (|2.19[> it is easy to compute 
the optimal signal-to- noise ratio p for the signal (|2.16l) : 



p := \J (h\h) = ho\ 



T 



Sn{fcY 



(2.20) 



Making use of Eq. f|2.16[) and the identities 2 sin a cos a = Substituting Eqs. f|2.19[) and (|2.16p into Eq. (|2.5p . we 

sin 2a, cos 2 a = A (1 + cos 2a), sin 2 a = k (1 — cos 2a) we get the following formula for the likelihood ratio: 



lnA(x; hi, h 2 ,£) 



2 T 



hi (a: sin $(<;;£)) + h 2 (x cos $(£;£)) - ^{h\ + h\) 



Sn(fc) 

To maximize In A with respect to the parameters hi and h 2 we solve equations 

<91nA 



dhi 



= 0, i = 1,2. 



The unique solution to these equations reads 

hi — 2(xsin$), h 2 — 2(x cos $) 



(2.21) 



(2.22) 



(2.23) 



Replacing in Eq. (|2 . 2 1 1) the parameters hi and h 2 by their estimators hi and h 2 [given by Eqs. (|2.23l) ]. we obtain the 
reduced likelihood ratio which we call the J- '-statistic: 



T(x; g) := In A(x; hi, h 2 , £) = 



2T n 



S„(fo) 

It is rather easy to rewrite the ^-statistic in still another form 

ii+T /2 



x sin$(t;£)) 2 + (x cos$(t;£)) : 



(2.24) 



(fc)T 



J x(t) exp ( - iu>i ) J exp ( 



- i^o^r ) di 



ti-T /2 



(2.25) 



Thus the ^-statistic, up to a constant multiplication fac- 
tor, is the modulus squared of the Fourier transform of 
the product of the data stream x(t) by the exponential 
factor exp f — iu>i (t/T ) ) which depends on the spindown 
paramerer w\. 



For the construction of the bank of templates it is cru- 
cial to study the expectation value of the ^"-statistic 
(|2.24[) in the case when the data x contains some 
gravitational- wave signal h, i.e. 

x(t) = n(t) + h(t;6'), (2.26) 
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where 6 1 — {h[,h' 2 ,£') collects the parameters of the 
gravitational-wave signal present in the data [see Eq. 
(|2.16|) ], Let us denote this expectation value by Ei (the 
subscript ' 1 ' means here that the average is computed in 
the case when the data contains some gravitational- wave 
signal), so we have 

Er{.F(z; 0} = E{^(n(t) + h(t; 8')- £)}. (2.27) 

We want to obtain an approximate analytical formula 
for this quantity. Making use of the following approxi- 
mations 

(sin[$(t; i) + *(t; £')]) = 0, (cos[$(i; £) + $(t; £')]} = 0, 

(2.28) 

after some computation we obtain 

E 1 {F(x; £)} = 1 + \ P 2 (< sin [*(t; - $(t; £')] >" 

+ <cos [$(i;0-*(t;O]) 2 ), 

(2.29) 

where p is the signal-to- noise ratio from Eq. (|2.20l) . The 
right-hand side of the above equation can be rewritten in 
terms of the autocovariance function Co of the ^-statistic 
(computed in the case when the data contains only noise); 
it is defined as 

Co(€,0 :=E{[jr(n;0-m (C)][-F(n;0-TOo(0]}, 

(2.30) 

I 



where mo is the signal-free average of J-: 

m (0:=E{7>;£)}- (2-31) 

In Sec. IV of Rcf. 10] it was shown that the autoco- 
variance function Co computed for the gravitational- wave 
signal of the form (|2.16p can be approximated by 

+ ( cos [<!>(*; £)-$(*;0]> 2 , (2-32) 

therefore the expectation value (|2.29p can shortly be writ- 
ten as 

E 1 {J-(x;O} = l + ^ 2 C (^|')- (2-33) 

The phase <3> [see Eq. (|2.1ip ] of the gravitational- wave sig- 
nal (|2.16[) depends linearly on the parameters £, therefore 
the autocovariance function (|2.32j) depends only on the 
differences between the parameters £ and 

Cfa&O = (sin$(f;^-^)> 2 + (cos<I>(t^-r)) 2 . (2.34) 

If one introduces t := £ — one can thus write 

C (t) S (cos$(t;r)) 2 + (sin *(t; r)) 2 . (2.35) 

Let us note that the function Co attains its maximal value 
equal to 1 for r = 0. 

We will numerically compute the autocovariance func- 
tion (|2.35[) in two ways. First, the right-hand side of Eq. 
(|2.35p can be expressed (without any additional approx- 
imations) in terms of the Fresnel integrals: 



C e (r, Xi ) 



2\ui\ 



ujq + uji (2xi + 1) 
V2 7r|wi| 

ujq +oji (2xi + 1) 



-F„ 



/ o; + cji (gxj - 1) 



LUq +UJ 1 (gXj - 1) 

v/2 7r|wi| 



n 2 



(2.36) 



where we have introduced the dimensionless variable 

ti 

Xi ■= 77-T • 
^ o 

The Fresnel integrals are defined atQ 

■/2 



(2.37) 



F c (x) 



■(H 



dz. 



(2.38b) 



F s (x) 



sm 



/7T , 

\~1 



We can also compute the right-hand side of Eq. (|2.35p 
in an approximate way. To do this we expand (|2.35p in 
dz, (2.38a) Taylor series around r = up to terms quadratic in r. 

Making use of the obvious equalities 



Let us note that both F c and F s are odd functions: F c (-as) = . <9 2 $ 

-F c (x) and F s (-a;) = -F„(a:). $ (*! £ = °) = °> = 



6 



we get 

2 

C a (T, Xi) = 1 - J2 f (*)w Tk Tl > ( 2 ' 4 °) 

fc,i=l 

where T is the 2-dimensional reduced Fisher information 
matrix with elements equal to 



\ 9r fc Or/ 



<9$ 
<9r fe 



fc, Z = 1,2. (2.41) 



In terms of the dimensionless variable \i the Fisher ma- 
trix r equals 



f(xO 



12 6 
,6 Xi 180 3 , 



(2.42) 



In Fig. Q] we study the relation between the exact (|2.36[) 
and approximate (|2.40l) formulae for the autocovariance 
function. We have found (see the left panel of Fig. [T]) 
that the approximate formula underestimates the value of 
the autocovariance function. In the right panel of Fig. [T] 
we plot the isoheights of the approximate autocovariance 
function and the isoheights of the fractional difference 



AC := 



C e — C a 

a, 



100%. 



(2.43) 



One can see that always AC > and for C c > 0.75 the 
fractional difference AC < 4%. 

From Fig. Q] one can also see that whereas the iso- 
heights of the approximate autocovariance C a are per- 
fect ellipses, the isoheights of the exact autocovariance 
C c are closed curves of shapes very similar to that of 
ellipses. Therefore it is reasonable to study the value 
of the approximate autocovariance along the isoheight 
C = const of the exact autocovariance. We have done 
this for the several values of C . Along each C e = const 
curve the values of the approximate autocovariance are 
smaller than C e and they are almost the same. We have 
picked up the largest value out of them and these values 
are given in Table U together with the corresponding val- 
ues of C e and the fractional difference between C a and 
C . The following cubic fit, 



C e 



-0.823872 + 3. 15065 C a - 1.84055 C a 2 - 



0.514037 Cl, 
(2.44) 



I/. 44 

reproduces the relation C c = C c (C a ) with accuracy bet 
ter than 0.04% for 0.55 < C a < 1.00. 



III. BANKS OF THE TEMPLATES 



To search for the gravitational- wave signal (|2.14[) in de- 
tector's noise we need to construct a bank of templates 
in the space of the parameters {ojq,ujx) on which the T- 
statistic [given in Eq. (|2.24[) ] depends. The bank of tem- 
plates is defined by a discrete set of points, i.e. a grid 
in the parameter space chosen in such a way, that for 



TABLE I. The relation between the exact 
the approximate, Eq 



(|2T40 

function of the ^-statistic. 



Eq. UH3BJ, and 

formulae of the autocovariance 





C a 


( Cc " Ca ) / Cc 


0.55 


0.43767 


20.4% 


0.60 


0.51513 


14.2% 


0.65 


0.58769 


9.6% 


0.70 


0.65600 


6.3% 


0.75 


0.72057 


3.9% 


0.80 


0.78183 


2.3% 


0.85 


0.84012 


1.2% 


0.90 


0.89575 


0.5% 


0.95 


0.94897 


0.1% 



any possible signal there exists a grid point such that 
the expectation value of the .F-statistic computed for the 
parameters of this grid point is not less than a certain 
fixed minimal value. From Eq. (|2.33[) we see that this 
expectation value depends on the signal-to-noise ratio p 
and on the value of the noise autocovariance function Co 
computed for the intrinsic parameters £ and £' of the 
template and the gravitational-wave signal, respectively. 
In the rest of this paper we will approximate the autoco- 
variance function Co by means of the formula (12.401) . i.e. 
we will use the equality Co(£,£') = C a (£, £')■ 

The signal-to-noise ratio p we can not control, therefore 
to construct the bank of templates one needs to choose 
some minimum value C m i n of the autocovariance function 
C a and look for such a grid of points that for any point 
£' in the (£Jo,Wi) plane there exists a grid node £ such 
that the autocovariance C a computed for the parameters 
£ and £ is not less than C m - m . Because C a depends on 
£, only through the difference £ — £ , we require that 



c a (£ - O > c„ 



(3.1) 



By virtue of Eq. (12.401) this condition leads to the in- 
equality 



^ ffcj (6c 
fe,i=i 



•&)(6-#)<i-C a 



(3.2) 



which for the fixed £ is fulfilled by all points £' which 
belong to an ellipse with the center located at 

Wa want to find the optimal grid fulfilling the require- 
ment p.ip , i.e. the grid which consists of possibly small- 
est number of points. Thus the problem of finding the 
optimal grid is a kind of covering problem, i.e. the prob- 
lem to cover the (wo,wi) plane (or, in data analysis case, 
the bounded region of the plane) by the smallest num- 
ber of identical ellipses. The thorough exposition of the 
problem of covering n-dimensional Euclidean space by 
identical spheres is given in Chap. 2 of Ref. (27j . 

We restrict ourselves to grids which are lattices, i.e. to 
grids with nodes which are linear combinations (with in- 
teger coefficients) of two basis vectors. If the vectors (Po, 
Pi) are the basis vectors of a lattice, then a fundamental 
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FIG. 1. Left panel: isoheights of the autocovariance function computed by means of the exact formula (|2.36f) (black lines) and 
the approximate formula (|2.40p (grey lines); isoheights for the values 0.5, 0.75, and 0.8 are shown. Right panel: isoheights of 
the autocovariance function computed by means of the approximate formula (|2.40|l (grey lines) and isoheights of the fractional 
difference AC defined in Eq. (|2.43p (black lines). 
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parallelogram (for the case of more than two parameters 
it would be a parallelotope) is the parallelogram consist- 
ing of the points 

AqPo + AiPx, 0<Aq,Ai<1. (3.3) 

A fundamental parallelogram is an example of a funda- 
mental region for the lattice, which when repeated many 
times fills the plane with one lattice point in each copy. 
The quality of a covering can be expressed by the cover- 
ing thickness 6 which is defined as the average number 
of ellipses that contain a point in the plane. For lattice 
coverings their thickness can be computed as 

g area of one ellipse _ ^ 

area of fundamental region 

We assume that we are interested in such searches 
for almost monochromatic gravitational-wave signals for 
which the number of grid points in the parameter space is 
very large and the time needed to compute the ^-statistic 
for all grid nodes is long, so it is crucial to use in the 
computation the fast numerical algorithms. Because the 
computation of the .F-statistic involves calculation of the 
Fourier transform [see Eq. (|2.25|) ]. one would like to use 
the fast Fourier transform (FFT) algorithm. The FFT 
algorithm computes the values of the discrete Fourier 
transform (DFT) of a time series. The values of the DFT 
are defined for a certain set of discrete frequencies called 
the Fourier frequencies. It will thus be possible to use the 
FFT algorithm in computation of the ^-statistic, if the 
grid points will be arranged in such a way, that the fre- 
quency coordinates of these points will all coincide with 
the Fourier frequencies. We have constructed two dif- 
ferent families of grids which fulfill this requirement, de- 
tails of the construction are presented in Appendix [Bj In 
the construction we have explored observation that the 
shape of the autocovariance ellipse strongly depends on 
the value of the parameter Xi'- the larger \xi\ is, the more 
alongated (along the wo-axis) this ellipse is (see Fig. [TJ). 
But the usage of alongated enough ellipses suggests that 
the .F-statistic could be computed for smaller number of 
frequency values than the standard DFT algorithm com- 
putes. Therefore we have employed some modifications 
of the DFT algorithm. 

Le the data collected by a detector form a sequence of 
N samples 

x u , u=l,...,N, (3.5) 

and let the sampling period be At. Then the DFT algo- 
rithm calculates the Fourier transform of the data with 
the frequency resolution A/ = l/(NAt). The resolution 
of the dimensionless frequency parameter ujq [introduced 
in Eq. (l2J|) ] is thus 

Alo q = 2vrT A/ = 2tt, (3.6) 

because AT At = T a . It is possible to modify the DFT al- 
gorithm in such a way, that the frequency resolution (13. 6[) 



changes. In Appendix |A1 we consider two such modifica- 
tions: (i) zero-padding of the data, which makes the DFT 
more dense; (ii) folding of the data, which diminishes the 
frequency resolution. Therefore we study banks of tem- 
plates which are compatible with the frequency resolu- 
tions of the form 

Aw = 2 £ 7T, £ = 0,1,2,3. (3.7) 

As explained in Appendix [21 t = corresponds to zero- 
padding of the data (for N data points we add N zeros), 
£ = 1 is a pure DFT of TV-point data stream, for I = 2 
the data is folded two times and for £ = 3 the data is 
folded four times. 

Let (Po, Pi) be the basis vectors of a lattice we con- 
sider. As explained above we want to use the DFT algo- 
rithm, therefore we need a such bank of templates that 
all nodes can be arranged along straight lines parallel to 
the wo-axis. Moreover, the distance between neighboring 
nodes along these lines must be equal to the frequency 
resolution {37?]) of the DFT algorithm. To fulfill this con- 
straint we require that the vector Po has components 

Po = (Auj q = 2 e 7T, 0). (3.8) 

Let us also denote the components of the second basis 
vector Pi as 

Pi = (&j 0l tfwi). (3.9) 

We will call the grid orthogonal, if Po • Pi = 0, where dot 
denotes the usual Euclidean scalar product. The grid 
spanned by the vectors (|3.8|) and (|3.9I) is thus orthogonal 
if and only if Sloq = 0. 

In Appendix [B] we have constructed two different fam- 
ilies of grids fulfilling the constraint (|3.8p . Construction 
of the grids denoted by Gi^ (for £ = 0,1,2,3) is de- 
scribed in Appendix IB 11 the grids G2/ and G' 2 t (valid 
for £ = 1,2,3) are constructed in Appendix IB 2 al and 
the grids G2.0 and G' 2 are described in Appendix IB 2 bl 
The grids Gi^ and G' 2 e are orthogonal, whereas the grids 
G24 are non-orthogonal. 



IV. TIME PERFORMANCE OF THE GRIDS 

We have made a number of numerical simulations to 
study the performance of the grids constructed in Ap- 
pendix [Bj All computations were done for N = 2 19 = 
524288 data points. For the sampling period At = 0.5 s 
(i.e. for the Nyquist frequency equal to 1 Hz) this corre- 
sponds to around three days of data. The dimensionless 
first spindown parameter u)\ we have taken to be non- 
positive and in the range U\ <E (—3000; 0). The discrete- 
in-time version of the .F-statistic [its continuous-in-time 
form is given in Eq. (I2.25|) ] reads 
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F(z;0 = 



2Ai 2 



N 



/ i u-r 

exp -icji I Xi- - + — ^~ 



cxp — iwq- 



u- 1 



(4.1) 



where the discrete data stream points are defined as 
x u := x(ti - T /2 + (u- l)At) (u = 1, . . . , N). The 
DFT algorithm computes the sum present in Eq. (|4.1[) 
simultaneously for all dimensionless frequencies ojq from 
a discrete set. Let us denote these discrete Fourier fre- 
quencies by ujQ r , then (see Appendix A) 



w r = (r - 1) Acj , r = l,...,r max , (4.2) 



where Aujq is the frequency resolution of the DFT. Using 
in the parameter space a discrete grid spanned by the 
vectors (Po,Pi) [see Eqs. (|3.8[) and (|3.9p ] means that 
also the spindown parameter uj\ becomes discrete with 



possible values equal to 

u Xs = -(s- l)£u>i, s= 1, ...,s max . (4.3) 

In the case of non-orthogonal grids the point in the 
parameter space with coordinates (ujq = 0, uji s ) usually 
is not a grid node for s 0; the grid node with the 
smallest wo-coordinate is (mod(s 5u>o, Aojq), wis), where 
mod(m, n) is the remainder on division of m by n. It 
means that for the fixed value of s ^ 0, before computing 
the DFT one has to multiply the uth data point by an 
extra exponential factor exp[— imod(s Sljq, Aojq)u/N]. 

Collecting all this information together one can finally 
show that the F-statistic evaluated on the grid point de- 
fined by the integers (r, s) [introduced in Eqs. (14.2[) and 
(|4.3[) ] can be written as 



T(x; r, s) = 



2At 2 
(fo)T 



^ x u exp 



! 1 w — 1 

i (s - 1) 6u x [ Xi ~ g + 



— i mod(s Sojo, ^ w o)^y- 



x exp 



-i(r - 1) Aw, 



u-1 



^ !)**•) ^max; S — 1 



} * * * j "max* 



(4.4) 



In our numerical simulations we have used the FFT 
algorithn0 with different frequency resolutions given in 
Eq. (|3.7I) . In Table [XT] we have shown the results of sim- 
ulations in which we have studied, for different types of 
grids, the time needed to compute the F-statistic for all 
grid nodes and the number of Fourier transforms per- 
formed during the computationjfl We have also shown 
the thicknesses 9 of coverings for different grids. 

Let us remind that in two dimensions the optimal cov- 
ering (without any constraint imposed) is the hexagonal 
covering with thickness 6 equal to 2^/(3^3) = 1.2092. 
The grids presented in Table |TT] have all thicknesses in 
the range 1.2106 < 8 < 1.3990, so they are greater than 
the optimal hexagonal covering thickness by ~0.1% to 
~15.7%. From inspection of Table ILTl it is clear that the 
value of the thickness does not decide on efficiency of a 



2 We have employed the fftw 3.2.2 implementation of the FFT 
algorithm with the wisdom mechanism (within which the algo- 
rithm performs some tests to ensure the optimal FFT perfor- 
mance). 

3 Our codes were written in C and compiled with gec 4.3.4. We 
have used PC computer with Core 2 Quad 2.66 GHz processor 
and 4 GB of RAM. 



grid. We see that the most efficient grids are grids which 
involve the computation of the longest FFTs. Moreover, 
within the grids with the same value of I (i.e. with the 
same length of the FFTs), much more efficient are orthog- 
onal grids. The reason for this is the lack of the expo- 
nential term exp[— imod(s (Jwo, Awo)w/A] in the discrete 
form of the F-statistic, see Eq. (|4.4[) . 



The non-orthogonal grid GJj'q is identical to the grid 
G2,0) but the time-performances of these two grids are 
different, because of different way of handling, during 
the computation of the F-statistics, the exponential term 
cxp[— imod(s 5uj , Au)q)u/N] in Eq. (|4.4[) . The comput- 
ing time can be considerably reduced if the values of this 
term are computed in advance for all needed values of s 
and u and they are kept in RAM memory of computer 
(note that the exponential term depends on the quanti- 
ties Awo and Sloq defining the grid and does not depend 
on data). However, to apply this trick one has to reserve 
a large amount of RAM for the table with values of the 
exponential term. For I = (i.e. for the grid G' 2 ' ) this 
table requires about 2.5 GB of RAM. One can estimate 
that the table for t + 1 requires around 2 times more 
RAM than the table for I. 
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TABLE II. The results of numerical simulations showing the time performance of the grids constructed in Appendix|B] Besides 
the times needed to compute the ^-statistic [given in Eq. 1)4. 4|l ] for all grid nodes, the number of the FFTs calculated during 
the computation and the covering thicknesses 9 are also shown. All grids are constructed for C m i n = 0.75 (which corresponds to 
Cmin = 0.7737). Each time slot in the table is the arithmetic mean from 200 repetitions of the computation of the J-"-statistic 
for all grid nodes. 



1 




Gl,£ 








G24, G' 2il 






Time (s) 


No of FFTs 


6 




Time (s) 


r" 


No of FFTs 


9 





23.4 


318 


1.2287 


48.1 


23.1 


30.4 


315 


1.2186 


1 


31.3 


705 


1.3639 


81.7 


32.2 




723 


1.3990 


2 


38.2 


1273 


1.2325 


126.6 


38.8 




1276 


1.2351 


3 


63.0 


2676 


1.2954 


230.2 


59.1 




2501 


1.2106 



Appendix A: Modifications of the discrete Fourier 
transform 



The data from the detector form the sequence 

(x r ) = (xx,x 2 , ■ ■ ■ ,xjv), (Al) 

so N is the number of data points. The discrete Fourier 
transform (DFT) of the data (|A1[) is defined as 



N 

JC Q ^^^^^ <ZJ J* 

r=l 



s (r-l)(a-l) 



N 



2;ri- — -). s = l,...,N. 

(A2) 

The DFT defined above computes the Fourier transform 
of the data stream (|A1[) at frequencies 



f - 5-1 -of _ 1^ 
h ~ NAt ~ [ ' N : 



s = l,...,N, 



(A3) 



where At is the sampling period and /n is the Nyquist 
frequency. The frequency resolution of the DFT (| A2|) is 
thus 



1 



NAt 



Awn = 2tt. 



(A4) 



1. Zero padding 



Let us now consider the 2iV-point data stream (y r ) 
which consists of the original ./V-point data stream (x r ), 
Eq. (jAll) . supplemented by N zeros, 

{y r ) = (x 1 ,...,x N ,0,...,0). (A5) 

By virtue of formula (| A2|) . the DFT of the data (y r ) reads 

y s = x r exp f - 2?ri ^ - ^ - ^ ) , s = l,...,2N. 



r=l 



2N 



(A6) 

Making use of Eqs. (|A2I) and (| A3|) one easily sees, that 
the numbers y s can be interpreted as the values of the 
DFT of the original data stream (x r ), but computed now 
for frequences 

L = S ~} , s = l,...,2N, (A7) 
J 2NAt ' ' K ' 



so the frequency resolution of the DFT (|A6I) is 
1 



A/ 



Aluq 



(A8) 



2NAt 

Of course one can add more zeros to the data to obtain 
their DFT with the frequency resolution better than this 
given in Eq. (|A8I) . 

2. Folding of data 

Let us fold the iV-point data stream (x r ) to the (N/2)- 
point data stream (y r ), 

N 

y r := x r + x r+N/2 , r=l, (A9) 
The DFT of the data (y r ) reads, according to Eq. (|A2|I . 



N/2 
J2(x r 



Xr+N/2) exp I - 2-n-i 



: (r-l)(s-l) 



N/2 



1 N 

8 = 1,...,-. 



(A10) 



It is not difficult, employing the periodicity of the func- 
tion exp(i2) for real z, to rewrite formula (|A10[) in the 
form 



N 

Vs = ^ Xr eXP 
r=l 



2tti 



s 2(r-l)(a-l) 



N 



1 H 
2 • 

(AH) 



Again making use of Eqs. (|A2I) and (IA3|) one sees, that 
the numbers y s can be interpreted as the values of the 
DFT of the original data stream (x r ), but computed for 
frequences 



f. 



2(8-1) 
NAt 



N 

I,...,—, 



so the frequency resolution of the DFT (|A11[) is 
2 



NAt 



Awq = 4"7T. 



(A12) 



(A13) 



By folding the data p times (p — 1,2,.. .) one gets (N/p)- 
point data stream the DFT of which is the DFT of the 
original Appoint data but computed with the frequency 
resolution Aluq — 2 P x 2tt. 
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Appendix B: Construction of the grids in the 
parameter space 

1. Grids Gi 

In this appendix we construct a family of grids by con- 
struction of their fundamental regions of parallelogram 
shape. The fundamental parallelograms are always in- 
scribed into the ellipse of constant value of the autoco- 
variance function [this ellipse is given by Eq. (|3.2[) ] and 
their construction ensures that the constraint p. 711 is ful- 
filled. 



a. Construction of the fundamental parallelogram 

Let us denote the coordinates of the parallelogram ver- 
texes we are looking for by (u)q ,u)^), a — 1,...,4. 
Bases of the parallelogram we choose to be parallel to 



the ui\ axis, so the ujq coordinates of the vertexes can be 
written as 

4 a \k) = k~, a = 1,2, 4 a) (fc) = -fc|, a = 3, 4, 

(Bl) 

where fc > 0. Then the coordinates ui[ a ^ of the vertexes 
we obtain from the equations [here T denotes the trans- 
position of the row vector (u^ ,u>i)]: 

(J a \J 1 a) yf( x oi4 a \"[ a) ) T = i-cim> o = i,. ..,4, 

(B2) 

where r(xi) is the Fisher matrix from Eq. (|2.42[) and 
(the reason for introducing the extra subscript 'i', 
from 'initial', is explained in Appendix IB 1 bl below) is 
the minimum value of the autocovariance function. We 
can assume that C min £ (0; 1), then the solution of Eqs. 
dB2j) reads 



(1) _ -30 fc tt xi + VTE y/48 (1 - CV n ) (1 + 60 X f) - fc^ rm , 

(*,*)- 2(l + 60 Xl 2 ) ' ( ' 



„(2), v k s _ -30 fcTrxi - (1 - CVJ (1 + 60 X f) £gg 

Wl (Xi ' fe) ~ 2(l + 60 Xl 2 ) ' [ ' 

a J [ 3 \ Xi ,k) = -u[ 2) (Xi,k), J^(xi,k) = M 1] (Xi,k). (B3c) 



Area of the parallelogram with the vertexes given above 
is equal to 



c i t 1 ) ( 2 )^ i 



(B4) 



Making use of Eqs. (|B3|) we get the area S as a function 
of the parameters Xi an d k: 



(B5) 

We want to have the fundamental parallelogram of pos- 
sibly large area. Let us first fix the parameter Xi an d 
maximize the area with respect to the parameter fc. The 
only positive solution of the equation dS/dk — reads 



ko(Xi) = ~ ^6 (1-^X1 + 60*?). (B6) 

TT v 

We can compute the area of the parallelogram for this 
optimal value of fc: 

Smax = S( Xi ,fco(Xi)) = 24^15(1 - CVJ. (B7) 

It does not depend on the vaule of Xi- It can be shown 
that the value S m & x given above is also the global maxi- 
mum of the function S — S(xu k)- 



The dependence S = S(xi,k) is illustrated in Fig. [21 
It is easy to show that for any value of fc > fc m i n , where 



fc mi „:-- X /6(l-CV n ), 



(B8) 



we can find two values of Xi f°r which the maximal area 
Smax is achieved. These values one obtains solving the 
equation S(xu k) — S max . The solutions read 



.B^iE^gj, (B9) 

For fc = fcmin we get only one solution Xi(^min) = 0. 

b. Construction of the grids 

We want to use the FFT algorithm in the computa- 
tion of the ^-"-statistic for all grid nodes, therefore we 
need such a grid that (i) all grid points can be arranged 
along straight lines parallel to the wo-axis, and (ii) the 
distance between neigbouring points along these lines is 
fc?r, where fc = 2 e , I = 0,1,2,3 [see Eq. l|57|l]. It is 
not difficult to construct an orthogonal lattice fulfilling 
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4^10 




FIG. 2. In the left panel the area S of the fundamental parallelogram for the grid Gi is plotted as a function of the parameters 
Xi and k [see Eq. (|B5|) ]. In the right panel the part of the (xi, k) plane is displayed; here the shaded region is made of the points 
(Xi, k) for which the construction of the parallelogram is not possible, and the V-shaped solid line is made of these values of 
the parameters (xi, k) for which the maximal is achieved. 



these requirements. Let the lattice be spanned by the 
vectors (Po,Pi). The first basis vector P is parallel 
to the wo-axis and has components Po := (fc7r,0), so its 
length is equal to the height of the fundamental paral- 
lelogram. The second basis vector Pi is chosen to be 
parallel to the wi-axis and has length equal to the length 
of the parallelogram's base. The wi-component of the 
vector Pi can easily be obtained from Eqs. (|B3|) . One 
finds that Pi = (0,Slui), where 



Soji 



15^ /48 (i-cvjq + eoxn-fc 2 ^ 

1 + 60 X f 



(BIO) 



In the construction of the fundamental parallelogram 
presented above (in Appendix IB 1 aft we have ensured 
that the autocovariance between the center of the fun- 
damental parallelogram and any of its vertexes is equal 
to CJqJjj. In the bank of templates spanned by the vectors 
Po and Pi the minimum value [taken over all points in 
the (luq, wi)-plane] of the maximum (taken over all grid 
points) autocovariance between any point in the (ujq, wi)- 
plane and all the grid points is not less than CL^. It 
turns out that usually it is greater than C l min . Let us de- 
note this minimum value of the autocovariance by C m ; n , 
then C min ^ CLj n . 

To find the value of C m i n we should first find the 
Voronoi cell of the grid0 To do this let us first choose 
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FIG. 3. The relation, valid for the grids G\, between the 
initial minimum value C^in an d the final minimum value C m i n 
of the autocovariance. Plots for three different values of the 
parameter k are shown. The diagonal grey line is the plot of 
the relation C m m = C' min . 



4 With each point £ a of a lattice C we associate its Voronoi cell 
V(£ a ) which consists of those points of plane that are at least as 
close to £ a as to any other point £ b £ C: V(£ ) := {£ : |£ — £ tt | < 
|C - €bl for all 6}. 



two neighbouring grid nodes and find all points in the 
(cjo, wi)-plane such that the autocovariance between the 
point and one of the chosen nodes is equal to the auto- 
covariance between the point and the second node. All 
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such points belong to the border between neighbouring 
Voronoi cells. We repeat this construction for all nodes 
which are neighbours of the chosen node. Points where 
borders cross each other are vertexes of the Voronoi cell 
and at these vertexes the autocovariance takes its mini- 
mum value. In Fig. [3] we have shown, for different values 
of the parameter k, the relation between the initial min- 
imum value C^jj, and the final minimum value C m - ln of 
the autocovariance. 



2. Grids G 2 

Constructions of the grids Gi were inspired by the ex- 
istence of the optimal hexagonal covering by circles of 
the 2-dimcnsional Cartesian space K 2 and they can be 
treated as some deformations of this covering. To em- 
ploy the properties of the hexagonal covering we firstly 
translate, by means of a linear transformation, the prob- 
lem of covering the (wo, wi)-plane by identical ellipses to 
the problem of covering the (u)' Q , w^)-plane by unit circles. 

In the constructions of the grids described below we 
were guided by the two following features of the hexag- 
onal covering, (i) The fundamental region of the hexag- 
onal covering can be chosen to be a regular hexagon in- 
scribed into the circle; we demand that deformed cover- 
ings have fundamental region in the form of a polygon 
(usually nonregular hexagon) inscribed into the unit cir- 
cle, (ii) Let us denote by (P^P^) the basis vectors of 
the lattice covering in the (w , wJJ-plane (with the vec- 
tor P is parallel to the w -axis); we demand that the 
lattice points which lie along two neighbouring straight 
lines parallel to the vector are shifted with respect to 
each other by half of length of the vector P[ (this feature 
can be seen in Fig. @|. 

We start from constructing a linear transformation (de- 
scribed by a matrix M) which converts the ellipse of 
the autocovariance function into the circle of unit ra- 
dius. The matrix M transforms a point with coordinates 
(wo,wi) into the point with coordinates (ojg,w^): 

(w ,wi) T = M-(u a ,u 1 ) T . (Bll) 

Equation of the ellipse of the autocovariance function in 
the (w , wi)-plane reads [see Eq. (|B2|) ] 

(loo, tax) ■ f ■ (w , wi) T = 1 - C min . (B12) 

The linear transformation (|B11|) converts the ellipse 
(IB12[) into the circle of unit radius provided the matrix 
M fulfills the condition 

M T • M = l - — f. (B13) 

1 ^min 

The Fisher matrix T is symmetric and [what can eas- 
ily be shown by means of Eq. (|2.42l) ] it is strictly pos- 
itive definite, i.e. (wo,wi) ■ T ■ (wo,wi) T > for any 
(wo,o;i) 7^ (0,0). For such matrix T the equation (IB13|) 



can be interpreted as its Cholesky decomposition, which 
states that there exists the unique upperdiagonal ma- 
trix M fulfilling Eq. (IB13|) . In the rest of this subsection 
we will assume that the matrix M is the result of the 
Cholesky decomposition (so it is an upperdiagonal ma- 
trix). Let us also note that the matrix M depends on the 
parameters Xi an d C m - m . 

Let Po be the vector in the (loq, wi)-plane parallel to 
the wo axis and with length equal to kir, so its (cl>o,wi)- 
components are Po = (kn,0). After transformation to 
the (w , u>[)-p\ane this vector becomes P . Let us de- 
note its length by akn, thus its (uj' , w^J-components are 
(remember that the matrix M is upperdiagonal) 

P = (afc7r,0). (B14) 

The construction of the second basis vector P^ of the 
lattice in the (w , w^)-plane is described in details below. 
When this vector is found, we take its image in the trans- 
formation inverse to that from Eq. (IB11[) . it defines the 
second basis vector in the (ujq, wi)-plane, 

Pi := M^ 1 P[. (B15) 

For the lattices constructed below the components of 
the basis vectors Pi depend on the parameter X'n there- 
fore for different values of this parameter we get different 
lattices, but, as we checked, all these lattices have (for the 
fixed t) the same value of covering thickness. Moreover, 
we can always choose such value of Xi that the vectors P 
and Pi will be orthogonal. In Table ITVl we give the com- 
ponents of the basis vectors (Po,Pi) for non-orthogonal 
grids (for £ = 0, 1, 2, 3) defined by choosing Xi = 0. 
Table fVl contains the components of the vectors (Po, Pi) 
for orthogonal grids G' 2 t (for £ = 0,1, 2, 3; we have added 
here primes to the grid symbols to distinguish them from 
non-orthogonal grids of Tabic HVjl together with the val- 
ues of the parameter Xi chosen to make the basis vectors 
orthogonal. 

a. Grids G 2 ,i (valid for £ > 1) 

To find the second basis vector P[ in the (w ,wj)- 
plane, we make the following construction, which is illus- 
trated in Fig. 0J We plot a circle K with radius equal to 
half of the length of the vector P . The center of the cir- 
cle coincides with the center of the line segment spanned 
by the vector P . We inscribe a right-angled triangle G 
into the circle K, its hypotenuse is along the diameter of 
the circle (i.e. along the vector P ). The vertex O at the 
right angle of the triangle G has coordinates 



0= i-akw-^-iakTry-cji 2 ,^). (B16) 

The second basis vector P^ is chosen to be parallel to one 
of the sides of the triangle G, see Fig. [4] We demand that 
the ratio of the length L of this side and half of length 
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TABLE III. Basis vectors Po and Pi denning the orthogonal grids G\ t t [t = 0,1,2,3) for C m i n = 0.75 (which corresponds to 
Cmin — 0.7737). The covering thicknesses 9 of the grids G\ t e are shown in Table UU 



£ 


a 1 

^min 


Xi 


Po 


Pi 





0.68038480461 


0.06911969777 


(7T,0) 


(0, 9.4565868877) 


1 


0.71208037943 


0.28027281768 


(2tt,0) 


(0,4.2593984076) 


2 


0.68138053992 


0.57228732748 


(4tt,0) 


(0, 2.3567814138) 


3 


0.69684651733 


1.19594612580 


(ftr,0) 


(0,1.1211909237) 



TABLE IV. Basis vectors Po and Pi defining the non- 
orthogonal grids G2,e (& ~ 0, 1,2,3) for C m m = 0.75 (which 
corresponds to C^m = 0.7737). The grids are defined by 
choosing Xi = 0. The covering thicknesses 9 of the grids G 2: e 
are shown in Table ILTl 



I 


Xi 


Po 


Pi 








(7T.0) 


(1.57079632679490, 9.53468292515346) 


1 





(2tt,0) 


(1.91202057746303, 4.15260608129565) 


2 





(4tt,0) 


(2.65356333373361, 2.35185575858832) 


3 





(8jt,0) 


(2.92397504375649, 1.19974457106559) 



TABLE V. Basis vectors Po and Pi defining the orthogonal 
grids G' 2 t e (£ = 0, 1,2,3) for C m i n = 0.75 (which corresponds 
to C^ in = 0.7737). The values of the parameter \i chosen to 
make the basis vectors orthogonal are also given. The covering 
thicknesses 9 of the grids G' 2t are shown in Table fill 



t 


Xi 


Po 


Pi 





0.08237276158660 


(7T.0) 


(0, 9.53468292515345) 


1 


0.23021935382642 


(2tt, 0) 


(0, 4.15260608129565) 


2 


0.56414244879677 


(4tt,0) 


(0, 2.35185575858833) 


3 


1.21858231921793 


(8tt,0) 


(0, 1.19974457106559) 




-4 -2 2 4 



FIG. 4. Construction of the grids G 2t i (valid for I > 1) in the 
(w , wi)-plane. 



of the vector P[ is an odd positive integer, L/(\P[\/2) = 
2n + 1 {n = 1,2,...). For the fixed n this requirement 
fixes the both coordinates of the point O and thus it 
determines the vector uniquely. Usually we obtain 
several possible values of n. We choose this value which 
leads to lattice with the smallest covering thickness 9. 



b. Grid G 2 ,o (valid for I = 0) 

The construction of the grid G^o in the (wq, o;^)-planc 
is illustrated in Fig. [5l We start from constructing three 
circles of unit radii. All these circles have to cross each 
other at the same point. We build then the three vectors 
Vi, v 2 , and v 3 of unit length (|vi| = |v 2 | = |v 3 | = 1). 
The vector Vi begins at the centre of one of the circles 
and ends at the point which is common to all three circles; 
the vectors v 2 and v 3 are constructed in a similar way, 
see Fig. [5] The coordinates of the vectors vi , v 2 , and v 3 
are 

v x = (?,<?), v 2 = (p,- 5 ), v 3 = (0,l), (B17) 
where p and q are positive numbers fulfilling the condition 
p 2 + q 2 = l. (B18) 

Because (see Fig. [5]) 

P =vi+v 2 , (B19) 
by virtue of Eq. (|B14[) (taken for k = 1) we get 

p = ^an. (B20) 

Making use of Eqs. (|B18|) . (|B20[) and the equality (see 
Fig. [5]) 

P'i = vi+v 3 , (B21) 

one easily obtains the (co' , wi)-coordinates of the basis 
vector P' 1; 



Pi= {\™,l + sjl-^a) "J. (B22) 
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